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ABSTRACT 

One of the key observables of the reionization era is the distribution of neutral and ionized gas. 
Recently, Furlanetto, Zaldarriaga, & Hernquist developed a simple analytic model to describe 
the growth of HII regions during this era. Here, we examine some of the fundamental simpli- 
fying assumptions behind this model and generalise it in several important ways. The model 
predicts that the ionized regions attain a well-defined characteristic size R r that ranges from 
~ 1 Mpc in the early phases to > 10 Mpc in the late phases. We show that R c is determined 
primarily by the bias of the galaxies driving reionization; hence measurements of this scale 
constrain a fundamental property of the first galaxies. The variance around R c , on the other 
hand, is determined primarily by the underlying matter power spectrum. We then show that in- 
creasing the ionizing efficiency of massive galaxies shifts R c to significantly larger scales and 
decreases the importance of recombinations. These differences can be observed with forth- 
coming redshifted 21 cm surveys (increasing the brightness temperature fluctuations by up 
to a factor of two on large scales) and with measurements of small-scale anisotropics in the 
cosmic microwave background. Finally, we show that stochastic fluctuations in the galaxy 
population only broaden the bubble size distribution significantly if massive galaxies are re- 
sponsible for most of the ionizing photons. We argue that the key results of this model are 
robust to many of our uncertainties about the reionization process. 
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1 INTRODUCTION 

The reionization of hydrogen in the intergalactic medium (IGM) 
is the hallmark event of the first generation of luminous ob- 
jects, marking the time at which they affected every baryon in 
the Universe. In the past few years, astronomers have made enor- 
mous strides in understanding this tra nsition, but the em erging 
picture is complex. The rapid onset of Gunn & Petersonl 1 19651) 
absor ption in the spectra of z > 6 quasars l Becker et al. 2001; 
iFan et alfeoollwhite et all2003h indicates that reionization prob- 
ably ended at z ~ 6 (albeit with large variance among differ- 
ent lines of sight: Sokasian et al. 2003 [their Fig. 10]; Wyithe & 
Loeb 2004a; Oh & Furlanetto 2005). On the other hand, the de- 
tection of a large optical depth to electron scattering for cosmic 
microwave background (CMB) photons iKogut et all 2003) indi- 
cates that the process began at z > 15 (although only with rela- 
tively low confidence). Reionization thus appears to be a compli- 
cated process, a con clusion strengthened by several other s tudies 
fTheuns et alJl2002t IWvithe & Loebll2004bl: iMesinger & Haimanl 
l2004lMalhotra & Rhoadsl2004l) . 

For these reasons, untangling the puzzle of how the "twi- 
light zone" of reionization ended is one of the most exciting ar- 
eas of cosmology. A number of new techniques are being devel- 
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oped to study this epoch. Perhaps the most interesting focus not 
on when reionization occurred but on how it happened: the evo- 
lution of the distribution of ionized and neutral gas can teach us 
an enormous amount about how the ionizing sources interacted 
with the IGM. The most promising approach is to observe the red- 
shifted 21 cm transition of neutral hydrogen in the IGM, which 
would allow us to perform "tomography" of the entire reioniza- 
tion process, watchin g in detail as the HII regions grew and filled 
the universe (see iFurlanetto & Briggsll2004l for a recent review). 
If the technical challenges inherent to these observations (includ- 
ing terrestrial radio interference, ionospheric distortions, and the 
Galactic and extragalactic foregrounds) can be overcome, such ex- 
periments promise to reveal the entire history of reionization. A 
number of other techniques will complement these studies, includ- 
ing the kinetic Sunyaev-Zel'dovich (kSZ) effect from patchy reion- 
ization in CMB m aps lAghanim et alll996llGruzinov &HiJl998t 
iKnox et al J 1998t) . absorption spectra of high-reds hift quasars, and 
the sp atial distribution of Lya-selected galaxies ( Furlanet to et alJ 
12004 J) . 

The key observables in each of these techniques are the sizes 
of HII regions and their distribution throughout the universe. These 
bubbles presumably appear around the first protogalaxies; as more 
ionizing sources form the bubbles grow and merge, eventually fill- 
ing all of space. A consensus is now emerging that these bub- 
bles likely attain quite large ( > 10 comoving Mpc) scales through- 
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out most of reionization: analytic arguments suggest large-scale 
variations in the ionizing emissivity lBarkana"& Loeb 2004j|), and 
numerical simulations show that indeed the HII r egions r apidly 
grow to sizes compara b le to the simulatio n boxes {Gnedin l 200C ; 
i Razoumov et alj 120021 ICiardi et alJ 120031 : ISokasian et alTboO; , 



2004). IWvithe & Loebl ©(Ma) 

also argued that the ionized bub 



bles must a chieve size s > 10 physical Mpc at the end of reion- 
ization, and Cen 1 2005) showed that, if dwarf galaxies in present- 
day clusters were responsible for reionization, the HII regions must 
reach similar sizes. Measuring the sizes of these bubbles through- 
out reionization, as well as their locations relative to galaxies and 
quasars, can teach us a great deal about structure formation in the 

early universe. 

Recently, iFurlanetto et alj J2004cl hereafter FZH04) con- 
structed the first quantitative model for the distribution of 
bubble sizes throughout reionization. From such a model, a 
wide range of predic tions are possib le, including those for 21 
cm surveys tFjjrlangJtoetalJ [2004b), quasar absorption spec- 
tra jFurlanetto et alJ I2004J) . the kSZ effect IZahn et alJ I2005L 
| McOuinri^t^ni2bo5al) . and Lva galaxy survevs nFurhrneTtoeHil] 
2004a, 2005). The model is based on the excursion set formal- 
ism f amiliar from the dark matter halo mass function I Bo nd et alJ 
Il99ll) . It uses a simple photon counting argument to build HII re- 
gions around proto-clusters of galaxies. We review the basic model 
in Sj2] One key prediction for all the observables is that the bubbles 
attain a well-defined characteristic size. In i|2. II and 12.21 we study 
how this scale depends on the underlying source population and 
the cosmology. 

To construct an analytic model, FZH04 made a number of sim- 
plifying assumptions. These include spherical symmetry, a constant 
ionizing efficiency throughout the universe, a simple approach to 
the halo mass function, an approximate treatment of the inhomo- 
geneous IGM, and ignoring "shot noise" in the galaxy distribu- 
tion as well as the potentially bursty nature of star formation at 
high redshifts. While some of these factors (such as the IGM and 
spherical symmetry) are sufficiently complicated that they must be 
addressed numerically, the others are amenable to an analytic treat- 
ment. lFurlanetto & OrJ <2005l) made an important step in this direc- 
tion by showing how recombinations affect the basic model: they 
have little effect until the bubbles reach their Stromgren radius rel- 
atively late in reionization. The main purpose of this paper is to 
examine how several other assumptions affect the major features of 
the model, especially those emphasised in £|2. 1 1 and to generalise 
it where possible. This is important because numerical simulations 
do not yet have the dynamic range to resolve the galaxies responsi- 
ble for ionizing the universe while simultaneously sampling a large 
enough volume to trace the growth of HII regions. These relatively 
simple analytic tests of the model offer the best hope of testing 
its applicability, although complete simulations will eventually be 
necessary to understand the entire process. At the same time, we 
will take advantage of the greatest strength of analytic models - 
their flexibility - to examine how basic properties of the galaxy 
population affect reionization and hence to isolate those properties 
to which observations will be most sensitive. We will find that the 
FZH04 model is surprisingly robust: its key features remain valid 
over a wide range in parameter space. Analytic models will thus 
play an important role in interpreting the evolving topology of neu- 
tral and ionized gas throughout this epoch. 

We consider several modifications to the basic formalism in 
fj3| We first allow a different halo mass function and a mass- 
dependent ionizing efficiency. We show that, although such mod- 
ifications can introduce measurable differences to the bubble distri- 



bution, they do not affect the qualitative aspects of the model. Then, 
in il3.3l we investigate whether stochastic fluctuations in the galaxy 
population can significantly affect the bubble size distribution. In 
fj4|we show how these factors affect the bubbles if recombinations 
are included. We give some examples of how our results affect ob- 
servables in J3] and we conclude in j6| 

In our numerical calculations, we assume a cosmology with 
Q m = 0.3, n A = 0.7, Q b = 0.046, H = lOO/ikms -1 Mpc" 1 
(with h — 0.7), n = 1, and as = 0-9, consistent with the most re- 
cent measurements lSpe rgeletail2 003). We use the transfer func- 
tion of Eisenstein & Hu 1 1998) to compute the matter power spec- 
trum. Unless otherwise specified, we quote all distances in comov- 
ing units. 



2 THE FZH04 MODEL 

The basic motivation of the model is to describe how HII regions 
grow and overlap around overdensities containing galaxies. We be- 
gin by assuming that a galaxy of mass m ga i can ionize a mass 
C^gai, where ( is a constant that depends on (among other things) 
the efficiency of ionizing photon production, the escape fraction of 
these photons from the host galaxy, the star formation efficiency, 
and the mean number of recombinations. Each of these quantities 
is highly uncertain (e.g., massive Populati on III star s can dramat- 
ically increase the ionizing efficiency: iBromm et aljEoOll) . but at 
least to a rough approximation they can be collapsed into this sin- 
gle efficiency factor. The criterion for a region to be ionized by the 
galaxies contained inside it is then / co n > where / co n is the 
fraction of mass bound in haloes above some m m i n . Unless other- 
wise specified, we will assume that m m i n = r/14, corresponding to 
a virial temperature T v i r = 10 4 K, at which point hydrogen line 
cooling becomes efficient. In the extended Press-Schechter model 
lLacev & C ole 1993), this places a condition on the mean overden- 
sity within a region of mass m, 

5 m > S x (m, z) = S c {z) - V2K ( [al in - a 2 {m)] 1/2 , (1) 

where Kq = erf _1 (l — C _1 ), u 2 {m) is the variance of density 
fluctuations on the scale m, cr^m = a 2 (m m i n ), and 5 c (z) is the 
critical density for collapse. The global ionized fraction is Xi = 
C/coii, g , where / co li,g is the mean collapse fraction throughout the 
universe. 

FZH04 showed how to construct the mass function of HII 



region s from 5 X in an analogous way to the halo mass function 
( Bond et aljl99ll) . The barrier in equation Q is well approximated 
by a linear function in a 2 , 8 X « B(m, z) = Bo + Bia 2 (m). In 
that case the mass function of ionized bubbles (i.e., the comoving 
number density of HII re gions with masses in the range w± dm/2) 
has an analytic solution I Sheth 1998; McOuinn et al. 2005a): 



nt(m) dm 



dlna 



x exp 



d In m 

B 2 (m,z 



2a 2 (m) 



B (z) 
a(m) 

dm, 



(2) 



where p is the mean density of the universe. This should be com- 
pared t o the analogous mass fu nction for virialised dark matter 
haloes <Press & SchechteJ 19741) : 



nh(m) dm 



d In a 



dlnm 



S c (z 



cr(m) 
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x exp 



2a 2 (r 



dm. 



(3) 



Note the close structural similarity between these two functions. 
Note also that, while every dark matter particle must be part of a 
virialized halo, only a fraction Xi should be incorporated into bub- 
bles. The integral of equation J2j is slightly different from this ex- 
p ected value because of t he linear barrier approximation (see eq. 1 1 
in McOuin n et all2005al) . but the error is typically < 1%. 

Several key points of the FZH04 model deserve emphasis; for 
concrete examples of these properties see, e.g., Figure [2] below. 
First, the bubbles are large: tens of comoving Mpc for Xi > 0.65. 
Furthermore, the bubbles attain a well-defined characteristic size 
R c at any point during reionization. The distribution of bubble sizes 
also narrows considerably throughout reionization. Finally, nb(m) 
varies only weakly with £ for a fixed Xi (or equivalently it is nearly 
independent of the redshift of reionization). In the remainder of this 
section, we will examine the physical origin of these properties. 



2.1 The characteristic bubble size 

The characteristic bubble size R c plays an analogous role to the 
"break" mass in the halo mass function. This is simply the 
peak of m 2 rih(rn), corresponding to the halo size that contains the 
most fractional mass. If we assume that |d In a/d In m\ ^constant 
[which would follow from a pure power-law matter power spectrum 
P(k) oc k n ], we find 



5 c (z). 



(4) 



In other words, the characteristic halo mass is simply the point at 
which the exponential cutoff sets in. Physically, it is the scale at 
which a "typical" 1-er density fluctuation will collapse on its own; 
this scale grows with time because S c (z) decreases. 

We can define R c in exactly the same way. Again assuming a 
power-law matter power spectrum, we have 



a 2 (R c ) 



-! + (! + AB qB( 
2B 2 



?R?W 2 



(5) 



The physical meaning is somewhat obscured in this case, but it has 



two simple limiting forms. First, if 4B B 1 <C 1, a(R c 



B . 



This occurs when Xi ~ 1 (Bo must be small if most of the universe 
is to lie inside HII regions). In this regime, the characteristic scale 
is simply the mass at which a typical \-a trajectory crosses the 
barrier. This is completely analogous to ro*. The final step, then, is 
to understand the meaning of Bo. From equation Q, it corresponds 
to the self-ionization condition at m = oo: 



C/coll(<5 = Bo,cr = 0) =1. 



(6) 



In other words, it is the overdensity required for the sources inside 
a given large region to ionize that region. Thus, R c is simply the 
scale on which a typical l-a density fluctuation contains enough 
galaxies to ionize itself. 

In the opposite regime, when 4Bq B\ ^> 1 (early in reioniza- 
tion), we have a(R c ) ~ Bo /Bi . In this case, the simple analogy to 
m-t, is less helpful, because R c is determined from the set of atypi- 
cal trajectories that actually cross the barrier (constituting a fraction 
Xi of all possible trajectories). The increasing barrier selects a pre- 
ferred scale from this set that depends on both Bo and the slope of 
the barrier. 

Given Xi and z, we can invert equations {6) and to find R c . 
Figure^ compares the result (solid curves) to the exact R c (dashed 
curves) at z = 14, 10, and 8 (from top to bottom). The solid curves 
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Figure 1. (a): Characteristic bubble size. Dashed curves show the exact R c 
for z = 14, 10, and 8 from top to bottom. Solid curves show the estimate 
of eq. |3J. (b): Ratio of the estimated characteristic size to the true value. 
The solid curves use eq. while the dotted curves include only linear bias 
corrections feq.l9l. 



in Figure^ show their ratio. Our simple estimate matches the true 
result quite closely, to within a factor of two for Xi > 0.2. (The cut- 
off at small x% occurs because the characteristic mass approaches 
C^min, the smallest physically allowed HII region.) The estimate 
is not perfect because we have neglected the curvature of the cold 
dark matter (CDM) power spectrum; it worsens at smaller Xi be- 
cause the trajectories sample a larger range of masses and also at 
lower redshifts because the galaxy distribution is more evolved. 

We can go one step further in order to understand how R c 
depends on the ionizing source population. To linear order we can 
write 



/coil (5, f = 0) = / Ammn h (m) [1 + S z b L (m)], 



(7) 



where S z is the physical overdensity at redshift z and 6l (m) is the 
Lagrangian linear bias l Mo & White 1996): 

S 2 (z)/a 2 (m) - 1 



b L (m) 



5° c 



with <5° = S c (z — 0). Then equation {5J becomes 



B, 



(8) 



(9) 



D(z) b cS ' 

where D(z) is the linear growth factor [D oc (1 + at high 
redshifts] and 



dmml)i(m) rih(m). 



(10) 



P/coll,g 

The dotted curves in Figure 0> show the results of this approx- 
imation (scaled to the exact solution). The poor accuracy of the 
linear theory estimate may seem surprising given that / C oii, g -C 1. 
One reason is that finite size effects suppress the galaxy population 
when R c < 10 Mpc. But a more important reason is cosmological: 
the shape of the CDM power spectrum is such that a 2 oc i£ - ( 3+n ) 
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(approximately), with n — > 1 for R — > oo and 71 — > — 3 for i? — > 0. 
The bubbles sample a wide range of effective slopes as reionization 
proceeds; we have n ~ —2.3, —1.8, and —0.6 at R c = 1, 10, and 
100 Mpc. Thus, early in reionization the radius is extremely sen- 
sitive to a; in this regime, small errors in Bo correspond to large 
errors in R c . At the later stages, however, R c becomes more robust 
to the approximations and the estimate is more accurate. 

Although the b e g approximation does not provide a fully satis- 
factory estimate for R c , it is accurate to about a factor of two when 
Xi > 0.5. It also provides an important clue to the mechanism fix- 
ing the characteristic scale: for a given Xi, the crucial parameter is 
D(z)b c ff. More biased galaxy populations will drive reionization 
toward larger scales because they cluster more strongly. The galax- 
ies give a substantial "boost" to each mass fluctuation and allow 
larger scales to cross the reionization barrier. Thus, weighting mas- 
sive galaxies more heavily will generically create larger bubbles. 
Furthermore, equation {5J explains why n&(m ) is only w eakly de- 
pendent on redshift for a fixed Xi : as noted by lOhlj 19991) . the halo 
bias is such that the galaxy correlation length evolves fairly slowly 
from z ~ 6-20. A similar factor D(z)b e g determines R c here. 

Finally, we note that the mere existence of R c does not explain 
why rib(m) has a small-scale cutoff (the halo mass function, for 
example, becomes a power law for m <C rru). In the excursion set 
formalism, the cutoff occurs because 8 X increases linearly with a 2 , 
so that random walks are not guaranteed to cross the barrier at any 
scale. Physically, the barrier increases because the collapse fraction 
is suppressed as we move to smaller scales, because fewer modes 
are available and because of finite size effects. 



2.2 The width of the bubble distribution 

The other characteristic we wish to explain is why rib(m) narrows 
as Xi — > 1. To understand this, we first consider nt(m) if P(k) oc 
k n ; in this case Xi affects the logarithmic width of the distribution 
only through the barrier in the exponential factor (see eq.|2j- In fact, 
the distribution of bubble sizes narrows as Xi — > in such a cos- 
mology: when Xi is small, the dominant factor is exp[— B 2 / (2a 2 )], 
and the distribution decays exponentially at small masses. It actu- 
ally evolves somewhat more slowly than this single factor indicates, 
because the other terms in the exponential compensate: it even re- 
mains nearly invariant for Xi >0.5. We find that, for a power law 
P(k), the distribution decays exponentially with a characteristic 
non-dimensional length a 2 (R s ) / a 2 (R c ) = K > 3, where R s is a 
small-scale cutoff, with K increasing slowly as x% increases. This 
is the opposite of what we find in a CDM cosmology, where the 
distribution widens at early times. 

The fundamental reason for the narrowing at large R c must 
therefore be the shape of the CDM power spectrum. Again, we 
note that a 2 oc _R-( 3 +"> with n w -2.3, -1.8, and -0.6 at 
R c = 1, 10, and 100 Mpc. Thus, R s /R c ~ K~ 1/in+3) . When 
the bubbles are large, R s /R c w K~ 4 ; the two scales must lie 
close to one another because the variance changes rapidly with 
scale. On the other hand, when Xi <C 1 and R c ~ 1 Mpc, we 
find R s /R c « K~ 14 . Thus, R s <C R c , and the distribution is 
much broader at early times. 

Interestingly, this means that the underlying power spectrum 
has a substantial effect on the growth of HII regions during reion- 
ization. Given the large number of astrophysical uncertainties in 
any model of reionization, it seems unlikely that such behaviour 
can actually be used as a probe of the power spectrum. However, it 
does imply that the existence of a well-defined characteristic scale 



is a generic feature of any reionization model based on the under- 
lying galaxy fluctuations in a CDM model. 



3 GENERALISATIONS OF THE FZH04 FORMALISM 

We will now modify the FZH04 model in a number of useful ways. 
We will see that the features described in Sj2|are surprisingly robust 
to many of the simplifying assumptions. 



3.1 A mass-dependent efficiency C( m ) 

One obvious simplification of FZH04 was to take the ionizing ef- 
ficiency parameter £ to be independent of galaxy mass. In the lo- 
cal unive rse the star formation e fficiency /* oc m 2 ^ 3 in low-mass 
galaxies I Kauff mann et aj2003l) ; similar behavior could certainly 
occur for any of the parameters relevant to reionization, especially 
if they are driven by feedback. Here we will examine a general 
class of models in which £ = C( m ); we W1 H ta k e C K ra ° i m our 
numerical calculations. To do so we simply replace the condition 
C/coii = 1 with 



dm((m) mrih(m\8, M) 



(11) 



where nh{m\6, M) is the halo mass function within a large re- 
gion of mass M (or equivalently radius R) and smoothed over- 
de nsity 8. A ccording to the extended Press-Schechter model, this 
is lLacev & C old 19931) 



2 p 



n h (m\8,M) = W - 



x exp 



diner 



2 [Sc(z)-5] 



dlnm 



2 (m) 

-8} 2 



2[<7 2 (m) — a 2 (M)] 



CT 2( M )]3/2 

(12) 



The parallel with the original condition is more obvious if we define 
£ — f (5, M) to be the mass-weighted ionizing efficiency within 
the specified region; then equation i 1 1 i becomes C/coil = This 
condition allows us to find (numerically) the minimum 8 X for a re- 
gion to self-ionize, just as in equation Q, and hence the appropriate 
nb(m) (note that the barriers remain nearly linear in a 2 ). 

Figure|2|shows the resulting size distributions. The solid, long- 
dashed, and short-dashed curves have Xi = 0.16, 0.65, and 0.975. 
In panel (a), we show the distributions at z = 12. Within each 
set, the curves have a = 0, 1/3, 2/3, and 1, from left to right. 
The most important result is that R c increases with a. This is to 
be expected from equation |9j: massive galaxies are more highly 
clustered, so they drive reionization to larger scales. The degree 
of amplification decreases as Xi — > 1 because of the shape of the 
power spectrum, just as described in EI2.2I At early times, there is 
an additional contribution because the C( m ) excursion set barrier 
S x is steeper, which moves the peak to smaller a or large scales 
[see the discussion following equation The steepening occurs 
because more photons come from massive galaxies, whose abun- 
dances are more efficiently suppressed on smaller scales. Note, 
however, that R c never increases by more than a factor of a few: 
this is because the barriers tend to intersect fairly near the charac- 
teristic scale. Panel (bj shows the distributions at z — 9. 1 Clearly 
the differences are comparable to the higher redshift case. 



1 Note that we take m„ 
star-forming haloes. 



3.3r?i4 here in order to force £ > 1 in all 
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Figure 2. Bubble size distributions for ( cc m°. The solid, long-dashed, 
and short-dashed sets have Xf = 0.16, 0.65, and 0.975. (a): z = 12. 
Within each set, the curves assume a = 0, 1/3, 2/3, and 1, from left to 
right, (b): z = 9. We show a = and 2/3 in this case. 



Figure 3. Bubble size distributions for different halo mass functions at 
z = 8. From left to right, the sets of curves take Xi = 0.16, 0.4, 0.65, and 
0.95; within each pair, the solid and dashed curves use the Press-Schechter 
and Sheth-Tormen halo mass functions, respectively. 



3.2 The halo mass function 

While the Press-Schechter mass function is elegant and provides 
a reasonable description of halo populations in numerical simula- 
tion s, it is possible to impr ove the agreement at lower redshifts. 
The lSheth & TormeiJ 1 1999) mass function provides a better fit to 
the simulation results i Jenkins et all200ll) : 



71st(ti) dm 



.4 



Sc(z) 
x exp 



dlna 



d In m 



5l(z) 



1 + 



2a 2 (m) 



dm, 



(13) 



where A = 0.322, a = 0.707, and p = 0.3. Although originally 
motivated through fits to numerical simulations, this form also has 
some physical justification. Shet h et all 1200 ll) showed that the dy- 
namics of ellipsoidal collapse yield an excursion set barrier whose 
mass function is similar to equation 1 1 31 . At high redshifts, the 
mass function is less well-know n, but it seems to lie some where 
in between these two versions l Jang-Condell & Hernauist 2001; 



Reed et al. 2003). We will now generalise FZH04 to include the 
Sheth-Tormen mass function. 

The first effect of the new mass function is to change /coii, g : 
it is larger than the Press-Schechter value at these redshifts because 
the Sheth-Tormen function yields more high-mass objects on the 
exponential tail of the halo distribution. This alone does not, how- 
ever, imply that nt(m) will change: it simply means that the two 
mass functions require different mean efficiencies f. The size dis- 
tributions will change only if the local collapse fraction behaves 
differently; i.e. whether the solution to C/ co ii,st (S, cr) = 1 has 
different scale of density dependence. Thus, we need to compute 
the Sheth-Tormen mass function in an arbitrary region of over- 
density 8. Unfortunately, the form in equation 1 1 31 cannot easily 
be extended to finite regions; in contrast to flat and linear barri- 
ers, the relevant excursion set barrier changes shape after a shift of 



the origin. ISheth & Tormenl l2002t) found that one could approxi- 
mate the conditional mass function nsT(mj<5, M) through a Tay- 
lor series expansion of the barrier after a translation of the origin. 
We will use this approximate form; it should suffice to illustrate 
the sensitivity to uncert ainties in the mass function. 2 We note that 
Barkana & Loeb 1 2004) took the Sheth-Tormen abundances to vary 
with 5 in precisely the same way as the Press-Schechter function 
does; if this is a reasonable approximation (as it appears to be from 
their comparison to simulations), we would expect nt,(m) to be 
independent of the halo mass function. 

Figure [3] shows the resulting size distributions at z = 8. 
Clearly the halo mass function has only a subtle effect on ni,{m). 
The slight increase in R c occurs because the Sheth-Tormen mass 
function contains more massive, biased haloes. The corresponding 
excursion set barriers are nearly coincident; only when Xi is small 
do they differ appreciably. At higher redshifts, the differences are 
even smaller than shown in this figure (at these times the mass func- 
tions are so steep that slight differences in its shape cannot have any 
significant effect). Thus it does not appear that uncertainties in the 
halo mass function will significantly modify the bubble distribu- 
tion. 

For completeness, we also note that iMcOuinn et all (2005a) 
discussed some related issues. They found that the mass function 
- and even the filtering scheme - can have a significant impact on 
the duration of reionization because it affects / co ii,g (and its time 
derivative). Here we have shown that nt(m) is nearly independent 
of the halo mass function at a fixed Xi. 



2 We do note one subtlety with this method: the Taylor-series fit to the mass 
function provided by Sheth & Tormen 12002) does not agree exactly with 
equation 1131 and gives a different mean collapse fraction. We therefore 
renormalise all of the conditional / co ii,st by the ratio of the mean fractions. 
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3.3 Stochastic fluctuations in the galaxy distribution 

To this point, we have assumed that the mapping between the mean 
matter density and the galaxy density is exact, so that the ion- 
ized bubbles trace the large-scale structure in a deterministic fash- 
ion. However, in reality the galaxy distribution will fluctuate. This 
may be important, because the FZH04 model shows that - par- 
ticularly toward the end of reionization - slight overdensities in 
the galaxy distribution can strongly affect the bubble sizes. This 
is a ma nifestation of the "phase tran sition" nature of reio nization. 
ISheth & Lemsonl \ 19991) and ICasas-Miranda et alJ J2002I) studied 
this "stochastic bias" in some detail. They found that the variance 
of galaxy counts in simulations of the z = universe is nearly Pois- 
sonian in underdense regions, somewhat smaller in regions near the 
mean density, and somewhat larger in overdense regions. They ex- 
plained the variation through finite size effects and clustering. We 
will consequently assume that the halo number counts at each mass 
obey Poisson statistics. Most of the ionized regions have only mod- 
est overdensities, so comparison to the lower-redshift simulations 
suggests th at we may overestimate the variance by a factor of order 
unity; Casas-Miranda et alJ l2002ll found the discrepancy relative 
to Poisson to be no more than a factor of two. 

How will such fluctuations affect 7Zf,(m)? Physically, the 
FZH04 condition C/coii = 1 compares the number of ionizing pho- 
tons to the number of hydrogen atoms within a specified region. 
Thus the relevant figure of merit for estimating the accuracy of our 
prescription is the variance in the number of ionizing photons pro- 
duced inside a bubble, afj . Each galaxy ionizes a number of atoms 
equal to Ni, g — C,m g jmH, where m g is its mass and run is the 
proton mass. If all galaxies had the same mass, the variance would 
be a 2 N = (£m,g/mH) 2 fihV (assuming Poisson statistics), where 
fih is the mean galaxy number density and V is the volume of a 
region with mass M. Allowing for a range of galaxy masses, we 
have 



a 2 N (M) = V / dm( 2 (m) — n h (m\S x ,M). 



(14) 



Fluctuations in the galaxy distribution correspond to changes in the 
number of ionizing photons per baryon. We can therefore quantify 
their effects in terms of fractional fluctuations in the local ionizing 
efficiency (: 



S ( = a N /{Ni), 



(15) 



where (Ni) = M/rriH is the expectation value of the total number 
of ionizing photons. 

Figure [4] illustrates the importance of these fluctuations. For 
reference, panel (c) shows nt,(m) in the standard FZH04 model for 
several different Xi and redshifts. Panel (a) shows 8^ as a function 
of scale. For this fiducial model, the fluctuations are of order unity 
for 7? < 0.3 Mpc and fall below 10% for R > 3 Mpc. Because the 
bubbles pass this threshold fairly early in reionization, they should 
not significantly alter nb(m) unless Xi is relatively small. 

Interestingly, the fluctuations are nearly independent of both 
redshift and Xi. To understand this behaviour, it is useful to rewrite 
equation 1151 as 



1 



(16) 



VfihV (m) ' 

where we have assumed £ = constant and the averages are over 
the conditional mass function. The first factor is the variance 



cPoisson 



N. 



-1/2 



gal 



expected from pure Poisson fluctuations; the 
second describes the distribution of galaxies. Note that both fac- 
tors implicitly depend on the bubble mass M. Let us first focus on 




Figure A. (a): Relative fluctuations in the number of photons per baryon 8q. 

(b) : Ratio of fluctuations in f to those expected from number counts alone. 

(c) : Bubble size distributions according to the FZH04 model (neglecting 
stochastic fluctuations; these are included only for reference). Each curve 
takes f =constant. The solid, long-dashed, and short-dashed curves assume 
Xi = 0.4 at z = 8, 10, and 15, respectively. The long dot-dashed, short 
dot-dashed, and dotted curves show the same redshifts but for Xi = 0.8. 



5^ olsson and assume that all galaxies have the same mass. Then (for 
a given bubble size) N sal oc l/£ and S^° iBBOn oc ( 1/2 ; clearly in 
this simple case the fluctuations would evolve with redshift (and 

Xi). 

However, in Figure|4j> we show that, because of the mass func- 
tion, jP° lsson does not necessarily dominate unless the bubbles are 
small (so that iV ga i < 1). The range of galaxy masses increases with 
cosmic time, increasing (m 2 ) and hence 8^ relative to the Poisson 
expectation. Evidently, evolution in the shape of the mass function 
(specifically m m i n /?n*) compensates for the increased number of 
galaxies at lower redshifts, making 8^ more or less constant with 
redshift. 

Equation d 1 5 i also helps illuminate why 8q changes only 
slowly with Xi at a fixed redshift. If the bubbles always took the 
mean density, their interior mass functions would be independent 
of £ and both factors in equation d 1 5i would be independent of 
Xi. However, the bubbles actually sit at density S x , which is a 
function of £;; as a result the shape of the mass function in ion- 
ized regions changes slightly because bubbles are more atypical 
regions when Xi is small. This also explains why 8^/8^°' BBOn is 
constant in large bubbles: all have nearly the same S x ~ Bo and 
a <C 1, so nh(m\S x , M) approaches a constant shape. In this 
regime, S c oc iV g ~[ /2 oc R~ 3/2 . 

The fluctuations approach <5P olaaon i n regions where the mass 
function is steep; in that case most galaxies have m ~ m m i n . We 
see this in the z = 15 curves in Figure^. It also happens if we 
increase m m i n , pushing galaxies farther off on the exponential tail 
of the mass function. However, in that case 5q increases because 
iVg a i falls by an even larger factor. 

Another way to increase 5^ is to take £ oc m a , with a > 0. 
Figure|5|shows the fluctuations if a = 0, 1/3, and 2/3 (solid, dot- 
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R (comoving Mpc) 



Figure 5. (a): ionizing efficiency fluctuations 8^. (b): Bubble size distri- 
butions at z = 8 (ignoring stochastic fluctuations). The solid, dotted, and 
dashed curves assume f oc m a , with a = 0, 1/3, and 2/3, respectively. 
The left and right sets take %i = 0.4 and 0.8. 



ted, and dashed curves, respectively) at z — 8. Panel (a) shows that 
the fluctuations can be quite strong when the efficiency is a steep 
function of galaxy mass, exceeding 10% for J?<20 Mpc. In this 

case the second factor in equation I15i becomes (C^m 2 ) 1 ^ 2 / (Cm), 
which is large (~ 20 and 60 for a = 1/3 and 2/3) because it 
exaggerates the importance of massive galaxies. Evidently, so few 
sources contribute that our deterministic approximation may break 
down until quite late in the process. 

How will these fluctuations affect 7i(,(m)? Clearly our ana- 
lytic model breaks down if 8^ ~ 1, and numerical simulations 
will then be necessary to study the process in detail. Nevertheless, 
we can guess the qualitative effects. Consider a region that would 
nominally form a bubble of volume Vo (corresponding to a den- 
sity threshold S x ). First, suppose that it has fewer galaxies than ex- 
pected, with an efficiency f_ < £. In the excursion set formalism, 
this implies an effective barrier height S' x > 8 X . We could there- 
fore estimate the local bubble distribution by considering random 
walks beginning at (S x , Vo) and computing where they cross the 
8' x barrier. This is essentially equivalent to searching for the pro- 
genitor distribution at a slightly earlier redshift. iFurlanetto & Obi 
(2005) showed that most of the volume is contained in a single 
progenitor slightly smaller than Vo. Thus, negative fluctuations 
will smear the volume distribution by a factor Sv ~ 8(, where 
5 V = (V- Vo)/V . 

Positive fluctuations are slightly more subtle because bubbles 
can merge. Ignoring clustering, each HII region is surrounded by 
neutral gas with volume sa Vo/xi. Early in reionization, these thick 
walls prevent the bubbles from merging and 8v ~ • On the other 
hand, if 5^ > (1/xi — 1), the bubble will grow enough to touch 
its neighbor. For larger fluctuations, bubbles will grow both by ion- 
izing new gas and by consuming their neighbors; the volume fluc- 
tuation will then be 8v ~ 5(/(l — Xi). Because the extra sources 
need only ionize the thin walls between bubbles, a modest boost in 
the number of galaxies can greatly increase the bubble's size. Thus, 



near the end of reionization, small fluctuations can broaden the size 
distribution by significantly more than S( . Clustering will also in- 
crease the importance of merging. However, during this epoch the 
bubble sizes are most likely limited by recombinations. We now 
turn to that regime. 



4 BUBBLES AND RECOMBINATIONS 

FZH04 included recombinations in only the crudest possible way, 
incorporating the mean number of recombinations per hydrogen 
atom into £. This assumption of a spatiall y uniform recom bina- 
tion rate is a gross oversimplification l Miralda-Escude et al. 2000, 
hereafter MHR00): the IGM is itself inhomogeneous, so the recom- 
bination time varies across voids, filaments, and collapsed objects. 
MHR00 argued that (on a local level) reionization therefore pro- 
ceeds from low to high densities. 

IFurlanetto &~OM {2005 ) showed how to reconcile this picture 
with the FZH04 model in which reionization begins around large- 
scale overdensities (Ciardi et al. 2003; Sokasian et al. 2003, 2004). 
For a bubble to grow, the (local) mean free path of an ionizing pho- 
ton must exceed its radius. The mean free path is determined by the 
spatial extent and number density of dense, neutral blobs (where 
the recombination time is short). Thus, as a bubble grows, its radi- 
ation background must also increase to ionize each of these blobs 
more deeply. As a consequence, the recombination rate inside the 
bubble also increases; HII regions cease growing when this recom- 
bination rate exceeds the ionizing emissivity (see also MHR00). In 
other words, an ionized region must fulfill a second condition: 

_ A d/ co u(<5, R) . ,„ „ 
ebub = C ^ > A bub (H,d) (17) 

where ebub is the normalised emissivity and ^4bub(-R, <5) is the re- 
combination rate within the bubble (which can be evaluated for a 
given IGM model). The latter depends on the mean overdensity of 
the bubble through the average recombination rate and increases 
with the bubble radius. This is analogous to the FZH04 condition 

C fcoii > 1 except that it compares the instantaneous, rather than cu- 

I I i i 

mulative, number of photons. Furlanetto & Oh (2005) showed that, 

for reasonable models of the IGM density field, equation I17i halts 
bubble growth at a well-defined radius i? max ; any bubbles nomi- 
nally larger than this size fragment into regions with R ~ -Rmax- 
(Note, however, that this saturation radius describes the mean free 
path of ionizing photons, not the extent of contiguous ionized gas; 
obviously as ii ^ 1 ionized gas does fill the universe.) Unfor- 
tunately, calculating i? max requires knowledge of the IGM den- 
sity field during reionization, which is currently an open question. 
MHR00 present a numerical fit to the IGM density in cosmological 
simulations at z — 2-4 that can be extrapolated to higher redshifts; 
for that model, i? max > 20 Mpc (and much larger when z > 10). 
Thus, it only significantly affects nb(ra) when Xi > 0.8. However, 
their model explicitly smooths the IGM density field on the Jeans 
scale of ionized gas, which may not be appropriate during reioniza- 
tion. If smaller-scale structures such as minihaloes are common be- 
fore reio nization, recombinations may become im portant much ear- 
lier (e.g. jHaiman et alj200ll:IShapiro et alj2004liiiev et alj2005l) . 
For concreteness we will use the MHR00 IGM model in these ex- 
amples. 
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Figure 6. Comparison of the emissivity inside a bubble (solid curve) to its 
recombination rate (dotted curve) as a function of bubble size at z = 8 
(assuming Xi = 0.8). The shaded envelopes show the range of 5 e . The hor- 
izontally and vertically shaded regions take p = 0.01 and 0. 1, respectively. 



4.1 Stochastic fluctuations and recombinations 

The FZH04 photon counting argument uses the cumulative number 
of ionizing photons, so the relevant fluctuating quantity is the to- 
tal number of galaxies in a region. The recombination limit, on the 
other hand, depends on the number of active sources at any time. 
Even if the total number of galaxies is fixed, each could shine for 
only a fraction p of the Hubble time (we will assume that p is inde- 
pendent of mass for simplicity). The variance in the luminosity of 
each galaxy is then (1 — p)/p times its mean square value; fluctu- 
ations in the emissivity are 



P 



(18) 



In the limit of a small duty cycle, the emissivity fluctuations ap- 
proach the number count fluctuations divided by p 1 ^ 2 , which ac- 
counts for the reduced number of active sources. The factor of 
(1 — p) in the numerator forces the fluctuations to vanish as the 
duty cycle approaches unity, because in that case all galaxies are 
active. 3 

Fluctuations in ebub will affect the distribution of ionized and 
neutral gas within a bubble as well as the rate at which it grows. If 
St > 0, it can ionize dense blobs more deeply, but that makes little 
difference because most photons can already reach the edge any- 
way. If <5 e < 0, the dense blobs will begin to recombine and more 
photons will be trapped inside. However, the effects will not be se- 
rious unless ebub falls below ylbub; even in that case, the bubble 
only shrinks if the recombination time is shorter than the timescale 
of the fluctuations (i.e., in dense regions along the edge). Fluctu- 
ations would be most interesting near i? max , because they would 

3 We have ignored fluctuations in the total number of galaxies here, be- 
cause we wish to consider emissivity fluctuations in a given bubble whose 
size is determined by the true number of galaxies, not its expectation value. 
Including both sources requires the replacement (1— p)/p —* 1/p. 



Figure 7. Comparison of the emissivity (solid curves) inside a bubble to 
the recombination rate (dotted curve) as a function of bubble size at z = 8 
(assuming Xi = 0.8). The shaded envelopes show the range of <5 e . The 
bottom curve with vertical shading takes f =constant and m m ; n = 7-714. 
The middle curve with horizontal shading takes ( oc m 1 / 3 and m m i n = 
7714. The uppermost curve with diagonal shading takes m m - ln = IO7714 and 
f =constant. All the shaded regions assume p = 0.1. 



then introduce an intrinsic sc atter to the saturation radius, which is 
otherwise quite well-defined l Furlanet to & Ohl2005l) . 

Figure|6|illustrates these effects. The dotted curve shows Abub 
(in units of recombinations per baryon per unit redshift) as a func- 
tion of bubble scale at z — 8 assuming x% = 0.8. 4 The thick curve 
show the emissivity ebub (in units of ionizing photons per baryon 
per unit redshift). Their intersection gives ii max . The shaded en- 
velopes delineate £bub(l ± 8 e ), illustrating the range of emissiv- 
ity fluctuations. The horizontal and vertical hatched regions take 
p — 0.01 and 0.1, respectively (note that the age of the universe is 
tu ~ 6.3 x 10 8 yr at z — 8). If the duty cycle is set by the dynam- 
ical time of galaxies (e.g., because active phases follow mergers), 
P ~ (p/Pvir) 1//2 ~ 0.075 and fluctuations are only important in 
small bubbles. On the other hand, if p< 0.01 (as may be appro- 
priate for quasars), fluctuations can be significant even in relatively 
large bubbles. We emphasise again, however, that the bubbles only 
shrink if both ebub < ^bub and the recombination time is short 
compared to ptn, which only happens for the mean density IGM 
when z > 10, even if p « 1. Also, the fluctuations remain relatively 
small near 7? m ax because of the enormous number of galaxies in 
these large bubbles. Thus, emissivity fluctuations should not intro- 
duce a larg e spread in the saturatio n radius, and many of the con- 
clusions of Furlanetto & Oh| 120051) . which relied on the existence 
of a sharply-defined characteristic radius, should remain valid. 
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4.2 Recombinations and the C( m ) barrier 

We now consider recombinations in the presence of a mass- 
dependent ionizing efficiency. This is easy to incorporate into 
Furla netto & Ohl 120051) : we simply replace the left-hand side of 
equation < 171 with d(£/ co n)/di. As mentioned in £13,11 weight- 
ing massive galaxies more heavily will increase the instantaneous 
emissivity at a fixed Xi because massive galaxies evolve more 
rapidly. Thus saturation from recombinations will thus occur at 
larger scales. FigureQillustrates this point. The bottom solid curve 
(surrounded by the vertical hatched envelope) takes £ ^constant, 
while the middle curve (surrounded by the horizontal hatching) 
takes £ oc m 1 / 3 . 

The ratio between the solid curves characterises the different 
rates at which the underlying galaxy populations evolve: the rele- 
vant haloes grow ~ 25% faster for £ oc m 1 ^ 3 . Thus the bubbles 
grow significantly larger when £ oc m 1//3 than when it is a con- 
stant parameter. For this example, i? max ~ 23, 63, and 180 Mpc 
for a = 0, 1/3, and 2/3. Of course, J? ma x also depends sensi- 
tively on the IGM density (we have used the MHR00 model here), 
and these quantitative predictio ns are thus not rel iable (causality 
will also limit the bubbles; IWvithe&Loebi r2004a). However, we 
have shown that whatever the underlying density structure, J? max 
increases significantly if massive galaxies drive reionization. The 
uppermost curve (surrounded by the diagonally-hashed envelope), 
which takes 772 m in = 10m4 and £ =constant, reinforces this point. 
Increasing the threshold mass by this factor makes the emissivity 
so large that 7? max > 300 Mpc (at least in the MHR00 model). 

Equation 1 1 81 shows that S e oc <5j, so the same factors that 
increase the number count variance will also increase that for the 
emissivity. The envelopes in Figure illustrate this: the relative 
fluctuations are largest for £ oc m 1//3 . They are nearly Poisson for 
m-min = 10m4, because most galaxies are close to the threshold 
and (m 2 ) 1//2 « m m i n . Note as well that by increasing ebut>, sce- 
narios relying on massive galaxies usually render emissivity fluctu- 
ations less important simply because tbub 2> ^4-bub on average. 



5 OBSERVABLE CONSEQUENCES 

The issues we have considered in this paper all affect the bubble 
size distribution during reionization. This is obviously one of the 
fundamental properties of any reionization scenario, and a number 
of ways have been proposed to study it. In this section we will 
briefly consider two of these techniques. For concreteness, we will 
contrast the original FZH04 model and a scenario with £ oc m 2//3 . 



5.1 21 cm measurements 

Several low-frequency radio arrays designed to detect high redshift 
(z ~ 6-20) fluctuations in the 21 cm emission from neutral hydro- 
gen are presently in the planning stages and/or under construction 
(e.g. PAST, LOFAR, MWA, and SKA). 5 These telescopes should 
answer many open questions about the reionization epoch. They 

4 The turnover at R < 1 Mpc occurs because we force each bubble to be 
ionized up to its mean density, even if the formal limit from the mean free 
path is smaller. ^^^^^ 

5 For more information, see |Pen et all 120051) . 
http://www.lofar.org/ http://web.haystack.mit.edu/arrays/MWA/ and 
http://www. skatelescope.org/ respectively. 



will measure fluctuations in Tb, the difference between the ob- 
served 21 cm brightness temperature and the CMB temperature T 7 
iFieldl 19591) : 

T * - 26 [(^)(T)] 1/2 
V Is J \ 0.022 J LVi2 m /i 2 / V 10 /. 

x (1 + 5) a; H i mK, (19) 

where T„ is the spin temperature and xm is the neutral fraction. 
Equation I19i ignores peculiar velocities, whic h tend to en hance 
the s i gnal for modes along the line-of-sight iBharadwai & Alii 
12004 iBarkana & LoeblEooa) . Modes perpendicular to the line- 
of-sight are, for the most part, unaltered by the velocity field. 
McO uinn et alj yOOSbj) show that when the HII regions are large, 
these redshift-space distortions are small on the large scales ac- 
cessible to upcoming experiments. We are interested in the 21 cm 
signal when the bubbles dominate, so we will ignore peculiar ve- 
locities. 

It is likely that X-rays from the first stars, shocks and black 
holes h eated the gas to temperat ures well above T~, early in reion- 
izatio n i5rl200lUVenkatesan et alfeOOlllChen & Miralda-Escudl 
120041) . At the same time, Lyman-continuum photons from the 
first stars coupled the spin temperature T s to the gas tem pera- 
ture th rough the Wouthuysen-Field effect 1 Wout huvserl 19521 : iField 
1958). While this process w as probably slower than the gas heat- 
ing, Ciardi & Madau (2003) argued that it becomes effective even 
when Xi < 1%. Therefore, we will take the limit T s 3> T 7 , such 
that the 21 cm brightness temperature at the location x is given by 

AT b (x)=T 6 (l-Xi[l + 6 x (x)])(l + S{x))-f b , (20) 

where Tb = fb /xm is the normalised mean brightness tempera- 
ture, xm = 1 — Xi is the global neutral fraction, 5 X is the over- 
density in the ionized fraction, and 5 is the matter overdensity (the 
baryons trace the dark matter on the relevant scales). 

The imaging capabilities of upcoming observatories will be 
limited, but statistical observations may stil l be quite powerful 
( Zaldarriaga et al. 20041 iBowman et"al]|2005l) . The simplest such 
statistic is the 3-D power spectrum (FZH04). Including terms up 
to second order in the perturbations 8, the brightness temperature 
power spectrum is 

PAT(k) = ft [x 2 H P SS + P XX - 2x„ P x5 + P x5x5 ] . (21) 

We construct P55 with the halo model ICoorav ^^hethj [2002]) 
and P x x and P x g with the methods described in McOuinn et al. 
(2005a). The fourth moment P x Sx5 follows from Pss, P xx , and 
P x s using the hierarchical ansatz. 

In FigureHl we plot [fc 3 PAT(fc)/27r 2 ] 1/2 at z = 12 for the 
original FZH04 model in which C, oc m° (solid curves) and for 
a modified model in which £ oc m 2//3 (dashed curves). For each 
model, we show Xi = 0.25, 0.4, and 0.7 (thick, medium, and thin 
curves, respectively). We see clear differences in the power spectra: 
both the peak (measuring R c ) and the amplitude differ between the 
models. The disparity decreases as a function of Xi, differing by 
about a factor of two for scales k < 1 Mpc at Xi = 0.25. This is to 
be expected from Figure|2|;: the differences between the bubble size 
distributions decrease with increasing ionized fraction, primarily 
because of the shape of the matter density fluctuations. While we 
have shown results for z = 12, the curves depend primarily on 
the bubble size distribution and show similar differences at other 
redshifts (see Fig.|2li). 

IBowman et alT J2005) find that the Mileura Wide- 
field Array (MWA) will be most sensitive to wavenumbers 
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Figure 8. 21 cm brightness temperature power spectra for x% = 
0.25, 0.4, and 0.7 (thick, medium, and thin curves, respectively). The 
solid and dashed curves take £ oc m° and £ <x m 2 / 3 , respectively. 

0.01 < k < 0.1 Mpc -1 (This is approximately true for most of the 
other planned telescopes as well.). These scales are larger than the 
effective bubble sizes for most Xi in our model. At such scales, 
equation <2U approaches 

Pat, (k) -> f b 2 [x% + b 2 - 2x H b] Pss , (22) 

where b is the volu me-averaged linear bias of the bubbles 
iMcOuinn et ai]l2005J) . In the (am' model, the bubbles are 
more highly biased. Intuitively, this is because they trace the most 
massive galaxies, which are of course also the most biased. Hence, 
for a fixed Xi, the large-scale signal increases in the £ oc m 2//3 
model. Conversely, in the small scale limit k > 1 Mpc, density 
fluctuations in the remaining neutral gas dominate the signal. It has 
the limiting form Pat,, (k) — -> T 6 2 x%x Pss ■ We indeed see that the 
curves with the same Xi converge at small scales; this regime con- 
tains little information about the ionizing sources. 

More generally, Figure [8] demonstrates how sensitive the 21 
cm signal is to the bubble size distribution. These observations will 
thus provide powerful constraints on ri{,(m). They should also be 
able to probe stochastic fluctuations in the galaxy population ( £I3.3> . 
which will decorrelate the Xi and the density fields. Other, more 
subtle properties (such as the halo mass function, °J3.2t may even 
be discernable in the distant future. However, we must note that the 
detailed power spectra will depend on other effects that we have 
not included in this simple analytic treatment - most important, 
we have assumed spherical bubbles throughout. These effects can 
be tested with numerical realisations of the models described here 
(such as those in IZahn et alj|2005l) or with full cosmological sim- 
ulations that also include radiative transfer; preliminary tests indi- 
cate that our analytic m odel roughly matches the numerical results 
fMcOuinn et all2005al) . except when Xi > .85 (at which point our 
picture of isolated bubbles breaks down). 

5.2 The kSZ effect 

The peculiar velocities of the HII regions imply that the scatter- 
ing of CMB photons off free electrons in the bubbles will impart 
a red- or blueshift to the sca ttered photons, creating hot and cold 
spots in the CMB lAghanim et aljll996l) . The power spectrum of 
temperature anisotropies produced by such scatterings, part of a 



broader class of anisotropies termed the kinetic Sunyaev Zel'dovich 
(kSZ) effect ISunvaev & Zerdovicblll98oll . have been calculated 
for several analytic models of reionization (e.g. iGr uzinov & Hj 
ll998tlKnox et alJll998tlSantos et alj|2003t IZahn et alj|2005t) . ReT 
cently, McOuinn et al. 1 2005a) used the FZH04 model to demon- 
strate that the kSZ anisotropies are quite sensitive to both the dura- 
tion of reionization and the bubble size distribution. 

In the original FZH04 model, l 2 Cf peaks at approximately 
the angular scale corresponding to R c when Xi ~ 0.5 (around 
/ = 3000) and has a substantial width in I owing to the distri- 
bution of bubble sizes. For I > 5000, the primordial component 
to the power spectrum becomes negligible owing to Silk damp- 
ing, providing a window in which we can observe these secondary 
anisotropies. The kSZ signal lies well below another secondary 
anisotropy, the thermal Sunyaev-Zel'dovich (tSZ) effect, at all an- 
gular scales. Fortunately, the tSZ anisotropies have a unique fre- 
quency dependence, which should allow them to be cleanly re- 
moved. The Atacama Cosmology Telescope (ACT) and the South 
Pole telescope (SPT) are designed to detect the kSZ anisotropies 
and should be operational as early as 2007. 6 We first summarise 
our fo rmalism for calculating the kSZ signal (see lMcOuinn et all 
2005a for more details) and then discuss how l 2 Cf changes if 
we allow £ oc m a . 

Thomson scattering of CMB photons off free electrons with a 
bulk peculiar velocity produces a temperature anisotropy along the 
line of sight n: 

AT f 

— (n) =o T \ dne~ T a(ri)n e {r],n) (n ■ v), (23) 

1 kSZ J 

where a is the scale factor, t(it) is the Thomson optical depth to 
the scatterer at conformal time r\, v(n, 77) is the peculiar velocity 
of the scatterer, and the electron number density is 

n e (77, n) = n e (rf) i, (77) [1 + 5 (77, n) + S x (77, &)]. (24) 

Here, n e is the average number density of electrons (free and 
bound). 7 

One might expect the dominant contribution to the kSZ angu- 
lar power spectrum to come at zeroth order in the perturbations 
{8, S x }, arising only from the velocity field. However, because 
only the component of the velocity along the line-of-sight con- 
tributes to the anisotropy and because it enters as an integral along 
the line-of-sight (eq. l23> . the signal averages over many troughs and 
peaks, making the net contribution effectively zero for the modes of 
interest. The lowest order contributions arise from modes perpen- 
dicular to the line-of-sight, whic h occur at s econd order in {5, S x }. 
iKaiserl 1 19921) and ljaffe & Kamionkowskl 1 19981) showed that at 
small scales and in the flat sky approximation the angular power 
spectrum of (AT/T)ksz can be written as 

C ' SZ = / % W( V ) 2 P^ sz (l/x, V ), (25) 

where x = yp — n, W = or n e (rjo) a -1 x 2 e~ T , and C Ma & Frvl 
l2002USantos et alJ2003l) 

P kSZ (k) = ^/^{(l^-k' 2 )^') 

6 For more information on these experiments, see 
http://www.hep.upenn.edu/act/act.html and http://astro.uchicago.edu/spt 
1 Note that eq. 1241 is not a formally rigorous perturbation expansion. 
However, it suffices for our purposes because the second order S X S term 
contributes a negligible fraction of the total kSZ signal. 
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x [P xx Qk - k'|) + 2P xS {\k - k'|)] } . (26) 

We have omitted the subdominant terms in equation <26> . and we 
use linear the ory to cal culate the power spectrum of velocity fluc- 
tuations P vv iMcOuinn et al. 2005a). 

Equation J26I also omits a Pss term that enters the sum within 
the brackets. This contribution is called the Ostriker-Vishniac 
(OV) effect; it has been studied in much more depth than the 
other "patchy" reionization terms (e.g., lOstriker & Vishniacll986t 
Jaffe & Kamionkowski 1998; Hu 2000), because it relies on much 
simpler physics. However, the other two terms dominate during 
the reionization epoch, so we do not include it in our analysis. On 
the other hand, the OV effect continues growing past reionization, 
and for reasonable reionization histories in the FZH04 model, the 
total OV signal tends to be ~ 1-5 times larger than the p atchy 
contributions around I = 5000 IZahn etaljEoollMcO uinn ct al. 
2005a). Fortunately, simulations should be able to calculate the 
OV anisotropies to reasonable precision, given a cosmology and 
a reionization history (they simply trace Pss for the baryons). In 
that case it can be isolated from the patchy reionization terms in 
which we are most interested. 

Figure [5] shows the kSZ angular power spectrum for models 
with ( oc m° (solid curve) and with £ oc m 2//3 (dashed curve). 
Both assume the same Xi(z) (see inset). In the £ oc m 2 / 3 model, 
the bubbles are larger throughout reionization. This is reflected in 
the power spectrum: it peaks at larger scales (I = 1000 rather than 
I — 2000) and is suppressed by ~ 30% relative to the constant £ 
model on small scales (/ > 3000). This occurs because fewer small 
bubbles appear. Unfortunately, the large scales on which the peak 
occurs will be washed out by the primordial anisotropies (the thin 
line in Fig.|9j. 8 In Figure|5| we plot the l-a ACT error bars (as- 
suming perfect removal of foregrounds as well as of the tSZ and OV 
anisotropies). If these contaminants can be adequately removed, the 
ACT should be able to distinguish these two models of reionization 
through the amplitude of the fluctuations on relatively small scales. 

We found that the bubble size distribution figures prominently 
in the angular distribution of power of the kSZ signal. However, 
the duration of reionization is primarily responsib le for de termining 
the amplitude of the kSZ signal l McO uinn et ai]|2 005a). The his- 
tory we have chosen is relatively narrow and has r = 0.095, which 
is about 2-a below the WMAP best fit value ofr = 0.17±0.04 
iKosut et all2003l) . The relative differences between the two mod- 
els should be similar regardless of the redshift of reionization or its 
duration, because they depend primarily on nt(m), which is nearly 
independent of redshift for a constant Xi. However, the amplitude 
of the signal could be up to a factor of five larger if reionization 
is extended, which would make these differences even easier to 
measure. The most important caveat is that the shape of the power 
spectrum on scales where the primordial anisotropies are negligi- 
ble changes relatively little between the £ oc m° and m' 2 ' 3 models. 
Distinguishing the size distribution from, for example, a slightly 
faster reionization history, may be difficult and will require con- 
sidera bly more d etailed modeling with simulations (c.f., Fig. 5 of 
iMcOuinn et all2005ii) . 



8 However, note that the kSZ effect on scales I ~ 1000-5000 will 
bias cosmological parameter measurements for future missions such as the 
Planck satellite. Thus the large scale behaviour is still detectable; it is just 
more difficult to isolate owing to additional degeneracies it suffers with 
the cosmological parameters that determine the primordial anisotropies 
ISantos et all2003tlZahn et all2005t) . 
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Figure 9. Angular power spectrum of the kSZ signal for models with 
£ oc m° (solid curve) and f oc m 2 / 3 (dashed curve). The inset shows 
the assumed ionization history. The thin solid line shows the primordial 
anisotropies, which dominate when I < 5000. The error bars show the 1-cr 
errors for the ACT experiment. 



6 DISCUSSION 

In this paper, we have taken a closer look at the key properties 
of the FZH04 analytic model for the growth of HII regions dur- 
ing reionization. First, we examined the origin of the characteristic 
scale R c . We showed that it has a simple physical interpretation as 
the scale on which a "typical" density fluctuation is large enough 
to ionize itself. The scale at which this condition is fulfilled de- 
pends on the bias of the underlying galaxy population, which gives 
a boost to the matter density fluctuations. Increasing the mean bias 
drives reionization to larger scales by allowing weaker density fluc- 
tuations (which appear on larger scales) to ionize themselves. 

Another key feature of the FZH04 model is that the bub- 
ble size distribution narrows throughout reionization while R c be- 
comes less and less sensitive to the input parameters. We showed 
that this primarily results from the shape of the underlying matter 
power spectrum: at Mpc scales, the effective power law slope is 
n ~ —2.3, while for R > 30 Mpc, n ~ — 1. In the former case, 
density fluctuations collapse at a similar time over a wide range 
of scales, making R c quite sensitive to the input parameters. For 
n ~ — 1, on the other hand, the scale dependence is much weaker. 

We then considered two modifications to the FZH04 formal- 
ism. One is the underlying halo mass function: FZH04 use the 
IPress & SchechteJ<1974l) distribution, which differs in detail from 
the mass f unction observe d in simulations of the low-redshift uni- 
verse |jenjdn^yii][2(XU]). We showed that the She th & Tormenl 
1 19991) parameterisation, which better matches simulations, barely 
alters the bubble mass function. The reason is that the two halo 
mass functions have similar dependence on the underlying density 
field. We also allowed the ionizing efficiency parameter £ to be a 
function of galaxy mass. This may be appropriate if, for example, 
feedback affected the star formation efficiency, IMF, or escape frac- 
tion inside galaxies. We showed that weighting massive galaxies 
more heavily increases R c because such galaxies are more highly 
biased. Under reasonable assumptions, this could increase R c by 
several times throughout the early and middle stages of reioniza- 
tion. 

This difference in the bubble sizes could be directly observ- 
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able in the near future. In particular, if massive galaxies drive reion- 
ization the 21 cm brightness temperature power spectrum increases 
by about a factor of two early in reionization on scales accessi- 
ble to the next generation of radio telescopes. The shape of the in- 
duced CMB angular power spectrum also changes at small scales, 
although it may be more difficult to extract information about the 
sources from this integrated measure because of degeneracies with 
the duration of reionization. Nevertheless, we have shown explic- 
itly that both of these techniques offer the exciting possibility of 
learning about the sources driving reionization - which may re- 
main hidden from us for many years - through measurements 
of the IGM on ~ 10 Mpc scales, which may be more accessi- 
ble. Another possibility is to constrain the HII distribution with 
the observed abundance and clustering of Lya-selected galaxies 
iFurlanetto et alJ2004iiEoo3) . 

Finally, we examined the conditions under which stochas- 
tic fluctuations in the galaxy distribution can significantly modify 
the bubble size distribution. Such variations are inevitable because 
galaxies do not trace the dark matter distribution precisely. Using 
a simple Poisson model (which is approximately correct at lower 
redshifts; Casas-Miranda et al. 2002), we showed that such fluctu- 
ations are typically < 10% once R c > 3 Mpc. In this case, they 
should have only a modest effect on the bubble sizes, except during 
the early stages of reionization. However, the fluctuations could be 
much more important if massive galaxies are responsible for reion- 
ization, because then the effective number density of sources falls 
dramatically. In such a scenario, numerical simulations of the full 
reionization process may be necessary to predict the bubble char- 
acteristics. 

One particularly interesting application of these results is the 
possibility of reionization by quasars. The bright quasars detected 
by the Sloan Digital Sky Survey do not seem able to reionize the 
universe at z ~ 6, prov ided that the luminosity function is not ex- 
tremely steep iFan et alj|2004l) . However, it remains possible that 
a population of relatively faint quasars could have made signif- 
icant contributions (though see iDiikstra et ai] Eo04). Our models 
suggest that such a picture would have a number of implications 
for the HII regions during reionization. First, if quasars preferen- 
tially traced massive galaxies, they would have produced charac- 
teristically large bubb les. This may be particularly relevant if, as 
recently suggested by iLidz et alj l20o"5l) . faint quasars (as well as 
bright ones) populated massive gala xies, or if quas ars grew most 
efficiently in massive haloes l Wvithe & Loeb 2002). Moreover, if 
quasars only populated massive galaxies, their stochastic fluctua- 
tions would be larger. This is particularly relevant because they 
also lik ely had short lifetimes with highly variable luminosities 
(Hopkins et al. 2005a b); Figure|6|shows that the resulting emissiv- 
ity fluctuations would cause the bubbles to grow in spurts. These 
global signatures o f quasar reionization complement the local sig- 
nature proposed by Za roubi & SiUdl2005h : the profile of ionization 
fronts could help distinguish hard and soft sources. 

We have found that the major properties of the FZH04 model 
of reionization are robust to several of the simplifications inherent 
to it. Most important, large-scale ( > 10 Mpc) inhomogeneities in 
the ionized fraction seem unavoidable. This is crucial to, e.g., 21 
cm measurements of reionization, because the next generation of 
telescopes will only be sensitive to structures several Mpc across. 
The existence of a characteristic scale is also generic to any model 
in which galaxies drive reionization. On the other hand, the FZH04 
model still contains a number of simplifications - crucially, the as- 
sumption of isolated spherical bubbles - that may be particularly 
problematic during the late stages of reionization. Cosmological 



simulations will still be needed to predict detailed observational 
signatures as well as to test the behaviour of our analytic models. 
Such tests are challenging because of the wide separation in scales 
between the galaxies and the HII regions. One promising hybrid 
approach is to implement the FZH04 model in numerical simu- 
lations or even simple gaussian random fields (Zahn et al. 2005). 
iMcOuinn et alj 12005 al) showed that the kSZ power spectrum from 
our fully analytic model and from numerical realisations match 
each other rather well, which bodes well for the utility of these 
simple analytic models. 

We thank M. Kamionkowski for helpful comments. This work 
was supported by NSF grants ACI 96-19019, AST 00-71019, AST 
02-06299, and AST 03-07690, and NASA ATP grants NAG5- 
12140, NAG5-13292, and NAG5- 13381. 
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